function out = inflVarianceRatio(outEstim,xSim)

model = outEstim.model;
setupFilter = outEstim.setupFilter;

% Solve for expected inflation up to max maturity+1
posInfl= find(strcmp(model.labely,'$\pi_t$'));
var0   = model.g0(posInfl,1);
varx   = model.gx(posInfl,:);
varxx  = model.gxx(posInfl,:,:);
varxxx = model.gxxx(posInfl,:,:,:);
maxMat = max(setupFilter.maturities);
[E_pai0,E_paix,E_paixx,E_paixxx]= projectionStepExpectedControl(var0,varx,varxx,varxxx,...
   maxMat+1,model,setupFilter.setupProj);

% Solve for expected policy rate up to max maturity
var0   = model.by0(1,:);
varx   = model.byx(1,:);
varxx  = model.byxx(1,:,:);
varxxx = model.byxxx(1,:,:,:);
[E_r0,E_rx,E_rxx,E_rxxx]= ...
    projectionStepExpectedControl(var0,varx,varxx,varxxx,...
    maxMat,model,setupFilter.setupProj);
% Adding r_cu to the solution
E_r0  = [var0  ;E_r0];
E_rx  = [varx  ;E_rx];
E_rxx = [varxx ;E_rxx];
E_rxxx= [varxxx;E_rxxx];

% One-step ahead expectations of nominal bond yields
nx       = model.nx;
E1_by0   = zeros(maxMat,1);
E1_byx   = zeros(maxMat,nx);
E1_byxx  = zeros(maxMat,nx,nx);
E1_byxxx = zeros(maxMat,nx,nx,nx);
for i=1:maxMat
    var0   = model.by0(i,:);
    varx   = model.byx(i,:);
    varxx  = model.byxx(i,:,:);
    varxxx = model.byxxx(i,:,:,:);
    [E1_by0(i,:),E1_byx(i,:),E1_byxx(i,:,:),E1_byxxx(i,:,:,:)]= ...
        projectionStepExpectedControl(var0,varx,varxx,varxxx,...
        1,model,setupFilter.setupProj);
end

% Decomposing news for the yield curve states in xSimusing simulated states
%xSim = outEstim.outFilter.outSim.xAll_cu;
T                = size(xSim,2);
nx1              = model.nx1;
E_infl_maxMat_cu = zeros(maxMat+1,T);
E_r_maxMat_cu    = zeros(maxMat+1,T);
yields_cu        = zeros(maxMat,T);
E1_yields_cu      = zeros(maxMat,T);
for t = 1:T
    if setupFilter.pruningOn == 1
        xf_cu = xSim(1:nx,t);
        if setupFilter.setupProj.appOrder == 1
            E_infl_maxMat_cu(:,t) = E_pai0 + E_paix*xf_cu;
            E_r_maxMat_cu(:,t)    = E_r0   + E_rx*xf_cu;
            yields_cu(:,t)        = model.by0 + model.byx*xf_cu;
            E1_yields_cu(:,t)     = E1_by0 + E1_byx*xf_cu;
        elseif setupFilter.setupProj.appOrder == 2
            xs_cu = zeros(nx,1);
            xs_cu(1:nx1,1) = xSim(nx+1:end,t);
            AA_cu          = kron(xf_cu,xf_cu);
            E_infl_maxMat_cu(:,t) = E_pai0 + E_paix*(xf_cu+xs_cu) + reshape(E_paixx,maxMat+1,nx^2)*AA_cu;
            E_r_maxMat_cu(:,t)    = E_r0   + E_rx*(xf_cu+xs_cu)   + reshape(E_rxx,maxMat+1,nx^2)*AA_cu;
            yields_cu(:,t)        = model.by0 + model.byx*(xf_cu+xs_cu) + model.Byxx*AA_cu;
            E1_yields_cu(:,t)     = E1_by0 + E1_byx*(xf_cu+xs_cu) + reshape(E1_byxx,maxMat,nx^2)*AA_cu;
        elseif setupFilter.setupProj.appOrder == 3
            xs_cu = zeros(nx,1);
            xs_cu(1:nx1,1) = xSim(nx+1:nx+nx1,t);
            xrd_cu = zeros(nx,1);
            xrd_cu(1:nx1,1) = xSim(nx+nx1+1:end,t);
            AA_cu           = kron(xf_cu,xf_cu);
            AAA_cu          = kron(AA_cu,xf_cu);
            BB_cu           = kron(xf_cu,xs_cu);
            E_infl_maxMat_cu(:,t)= E_pai0 + E_paix*(xf_cu+xs_cu + xrd_cu) + reshape(E_paixx,maxMat+1,nx^2)*(AA_cu + 2*BB_cu) ...
                + reshape(E_paixxx,maxMat+1,nx^3)*AAA_cu;
            E_r_maxMat_cu(:,t)    = E_r0   + E_rx*(xf_cu+xs_cu + xrd_cu) + reshape(E_rxx,maxMat+1,nx^2)*(AA_cu + 2*BB_cu) ...
                + reshape(E_rxxx,maxMat+1,nx^3)*AAA_cu;
            yields_cu(:,t)        = model.by0 + model.byx*(xf_cu+xs_cu + xrd_cu) + model.Byxx*(AA_cu + 2*BB_cu) ...
                + model.Byxxx*AAA_cu;
            E1_yields_cu(:,t)     = E1_by0 + E1_byx*(xf_cu+xs_cu + xrd_cu) + reshape(E1_byxx,maxMat,nx^2)*(AA_cu + 2*BB_cu) ...
                + reshape(E1_byxxx,maxMat,nx^3)*AAA_cu;
        end
    else
        x_cu   = xSim(1:nx,t);
        AA_cu  = kron(x_cu,x_cu);
        AAA_cu = kron(AA_cu,x_cu);
        E_infl_maxMat_cu(:,t)= E_pai0 + E_paix*x_cu + reshape(E_paixx,maxMat+1,nx^2)*AA_cu + reshape(E_paixxx,maxMat+1,nx^3)*AAA_cu;
        E_r_maxMat_cu(:,t)   = E_r0   + E_rx*x_cu   + reshape(E_rxx,maxMat+1,nx^2)*AA_cu   + reshape(E_rxxx,maxMat+1,nx^3)*AAA_cu;
        yields_cu(:,t)       = model.by0 + model.byx*x_cu + model.Byxx*AA_cu + model.Byxxx*AAA_cu;
        E1_yields_cu(:,t)    = E1_by0 + E1_byx*x_cu + reshape(E1_byxx,maxMat,nx^2)*AA_cu + reshape(E1_byxxx,maxMat,nx^3)*AAA_cu;
    end
end
E_realRates_cu      = E_r_maxMat_cu - E_infl_maxMat_cu;

%                     Et[infl_t+m] -                  - Et-1[infl_t+m+1]  
inflNews_per_period = (E_infl_maxMat_cu(1:maxMat,2:T) - E_infl_maxMat_cu(2:maxMat+1,1:T-1))*4;
inflNews            = repmat((1./(1:maxMat)'),1,T-1).*cumsum(inflNews_per_period,1);

realNews_per_period = (E_realRates_cu(1:maxMat,2:T) - E_realRates_cu(2:maxMat+1,1:T-1))*4;
realNews            = repmat((1./(1:maxMat)'),1,T-1).*cumsum(realNews_per_period,1);

yieldsNews          = (yields_cu(:,2:T)- E1_yields_cu(:,1:end-1))*4;
%yieldsNews         = (yields_cu(:,2:T)- yields_cu(:,1:T-1))*4;    %The martingale specification 
% Adding measurement errors to bond yields news
%yieldsNews          = yieldsNews+model.params.stdMr20*randn(size(yieldsNews));

premNews            = yieldsNews - inflNews - realNews;

% Output
if T < 1000
    out.realNews           = realNews;
    out.inflNews           = inflNews;
    out.premNews           = premNews;
    out.yieldsNews         = yieldsNews;
end

out.stdRealNewsPct     = std(realNews,0,2)*100;
out.stdInflNewsPct     = std(inflNews,0,2)*100;
out.stdPremNewsPct     = std(premNews,0,2)*100;
out.stdYieldsNewsPct   = std(yieldsNews,0,2)*100;


out.covRealInfl        = zeros(maxMat,1);
out.covRealPrem        = zeros(maxMat,1);
out.covInflPrem        = zeros(maxMat,1);
out.varReal            = out.stdRealNewsPct.^2;
out.varInfl            = out.stdInflNewsPct.^2;
out.varPerm            = out.stdPremNewsPct.^2;
out.varYield           = out.stdYieldsNewsPct.^2;
for i=1:maxMat
    tmp = cov([realNews(i,:)',inflNews(i,:)']);
    out.covRealInfl(i,1) = 2*tmp(1,2)*100^2;
    tmp = cov([realNews(i,:)',premNews(i,:)']);
    out.covRealPrem(i,1) = 2*tmp(1,2)*100^2;
    tmp = cov([inflNews(i,:)',premNews(i,:)']);
    out.covInflPrem(i,1) = 2*tmp(1,2)*100^2;
end
%out.varYieldCheck  = out.varReal + out.varInfl + out.varPerm + out.covRealInfl + out.covRealPrem + out.covInflPrem - out.varYield

out.RealVarRatio       = out.stdRealNewsPct.^2./out.stdYieldsNewsPct.^2;
out.InflVarRatio       = out.stdInflNewsPct.^2./out.stdYieldsNewsPct.^2;
out.PremVarRatio       = out.stdPremNewsPct.^2./out.stdYieldsNewsPct.^2;

matSelect    = [1 3 5 10]*4;
dataMoments  = [0.22 0.21 0.20 nan nan 0.46 0.44 0.46 nan nan 0.23 0.22 0.19 nan];
modelMoments = [ out.stdInflNewsPct(matSelect,1)' nan out.stdYieldsNewsPct(matSelect,1)' nan out.InflVarRatio(matSelect,1)'];
out.summaryTable   = [dataMoments' modelMoments'];

end

